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1.  Linear  non-autonomous  ODEs  with  parameters 


Example:  Sturm-Liouville  eigenvalue  problem 

y"(x)  +  u{x)  y(x)  +  A  y(x)  =  0 

plus  boundary  conditions,  or 


In  general:  need  to  solve  IVPs  of  the  form 


Y  =  A(t,\)Y,  Y(0)  =  Y0 

p 

A(t ,  A)  =  Ao(t)  +  A k  Ak(t) 

k= 1 


2.  Problem  features 


•  Need  to  solve  the  equation  for  many  different  parameter  values 

•  System  can  become  oscillatory  or  stiff 

•  Potential  is  slowly  varying,  and  typically  computed  numerically 

Idea 

•  Use  methods  that  work  well  for  near-autonomous,  oscillatory 
problems:  Magnus  integrators 

•  Collect  all  parameter  independent  parts,  and  evaluate  them  in  a 
precomputation  step 


3.  Neumann  series 


S  =  A(t)  S ,  5(0)  =  / 

can  be  written  in  integral  from  as 

(I  —  K)  o  S  =  S{t )  -  f  A(t)  5(r )  dr  =  I 

Jo 

Neumann  series  solution 

S(t)  =  (I  -  K)-1  o  I 

=  (I  +  K  +  K2  +  K3  +  . . . )  o  / 

rt  rt  nr\ 

=  I  +  A(t )  dr  +  /  A(ti)  /  A(t2)  dr2  dri  +  . 

Jo  Jo  Jo 

Convergence  criterion: 

|  A(t)  ||  dr  <  oo 


4.  Magnus  series 

Write  S(t)  =  ea^\  so  that 
a(t)  =  In  S(t) 

=  Ko/+(K2o/--(Ko/)2) 

+  (k3  o  7  -  t((K2  o  7)(K  O  7)  +  (K  o  /)( K2  o  /))  +  |(K  o  I)A 

+  •  •  • 

=  Si  +  S2  +  •  •  • 

where 


[A(n):  A(t2)]  dr2dri 


5.  Convergence  of  the  Magnus  series 


The  Magnus  expansion  converges  in  the  Euclidean  2-norm  provided 


f  ||y4(r)||  dr  <  — 

Jo  v 


where 


r2i r 

r0=  /  (2  + |r(l-cot(|r)))_1  dr  =  2.173737... 

Jo 

and  v  <  2  is  the  smallest  constant  such  that 

II [Aii ^2] II  <  Hl^ill  II ^2 II 

(Moan,  2002). 


Rough  estimate 


guarantees  convergence. 


f  ||7l(r)  ||  dr  <  1.09 
Jo 


6.  Toy  problem:  a  modified  Airy  equation 

y"(t)  +  (l  +  X)(t2  +  \)y(t)  =  0 


or 


with 


Ao(t)  — 


Y'  =  Ao(t)  +  \A1Y 
and  A\(t)  = 


0  1 

- 1 2  0 


'  0  1 
-1  0 


Note: 

•  For  A  large,  behavior  is  dominated  by  constant  frequency  fast 
oscillations 

•  Type  A  equation  according  to  Degani  and  Schiff  (2003) 

•  Quadratic  t  dependence  exposes  generic  leading  order  error  term 
for  moderate  A 


log  Error 


Global  error  for  A  =  1  at  t  =  10 


Global  error  as  a  function  of  A  at  t  =  10  with  N  =  213 


Global  error  as  a  function  of  A  at  t 


1  with  N  =  512 


Order  reduction  when  log  A 
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log1(J  Steps 


log  X  =  3.5 


log  A.  =  10 


7.  Classical  order 

Write 

A(t  -I-  /&)  =  cio  “I-  o\h  -I-  0,2 h2 

where,  for  the  modified  Airy  equation, 

(  0  1  +  A\  /  0  0\ 

a°  —  f—(i2  +  A)  0  )’  ai~{-2 to)’  “2 

Leading  local  truncation  error  for  forth  order  Magnus 


E™g 


h 5  fi 

^20  (2  ta°’  [a°’  a2l]  +  3  Iai>  iai’  a°]]  +  ta«’  [a«’  ta«’ 

4fi  i  \\  (2t(l \)(t2 -\- \)  (1  +  A) 

4U  +  Ajl  -bt2  +  A  -2t(l  +  \)(t2  +  \) 


720 

=  0{h5)0(  A3) 


+  e>(/ic 


8.  WKB  analysis 

On  the  interval  [t,t  +  h],  write 

A{t  +  h)  =  C{t)  +  H(t,h) 
where  H(t,  h )  =  0(h).  Assume  that 

C  =  XDX-1 

with  A  an  unitary  matrix,  i.e.  A^1  =  A*.  Factor  the  flow  map 

S(h)  =  X  eDhS(h) 

so  that 

S'(h)  =  A(t,  h)  S{h) 

5(0)  =  A-1 

with 

A(t,  h)  =  e~Dh  A-1  H  X  eDh 


9.  Modified  Magnus  -  RCMS 


The  rescaled  equation  can  be  used  as  a  basis  for  a  numerical  method, 
usually  in  the  form 

S\h)  =  e~ChHeChS{h) 

This  equation  can  be  solved  by  Magnus  expansion 

•  "Modified  Magnus  method"  (Iserles,  2002) 

•  "Right  correction  Magnus  series"  (Degani  and  Schiff,  2003) 

•  or  by  Neumann  expansion  (Iserles,  2004) 

This  approach  is  not  taken  here. 

•  Interferes  with  precomputation 

•  Implementation  issues 

•  Superior  behavior  only  for  very  large  A 


10.  Asymptotic  solution  for  modified  Airy 

For  the  modified  Airy  equation, 

(1)  A(t,  h )  =  exp(— Dh)  X a(t,  h )  X  exp (Dh) 

=  Ji(2 th  +  ft2)  (_  AT)  +  O  0)  . 

and  therefore,  the  first  term  of  the  Magnus  expansion  is 

h  =  l  A{t' hl)  dftl  =  F  (“2  +  j)  (o  -i)  +  °  (j) 


and  therefore 


S(h)  =  exp (Ch)  exp 


11.  Local  error  indicators  for  Magnus-4 


Recall  that 

where 


Smas{h)  =  eSl+S2 


rh 

si=  /  ( C  +  H{h1))dhi  =  Ch  +  Q{h 2) 

Jo 

82=  f  [  \c  +  H(h1),C  +  H(h2)]dh2dh1  =  0(h3\) 
Jo  Jo 


The  exponent  is  diagonalized  by  si  +  S2  =  Xm&gDmag(Xm&g)  1  where 

i\t2h5  fl  0 


Dmag  =  Dh  +  s!  - 


72  VO  -1 


+  h.o.t. 


Xmag  =X_ 


tti2  /o  0\ 

6^/2  V1  V 


+  h.o.t. 


12.  Local  error  indicators  II 


Altogether,  we  find  that 

Smag(h)  =  X  exp  (L>mag)  A-1 


tti2 

+  —  sin(/xma§) 
b 


1  O' 
0  -1 


+  h.o.t 


The  local  error  5(/i)  —  Smag(h)  has  leading  order  contributions 


E\oc  = 


4oc  = 


72 

fir 


th 


Xt2h 5  ^sin^™8  —  cos/imagN 
cos  /imag  sin  /imag 

1  M 

xo  -iy 

°  ^ 

1  0 


sin  (/ima§) 


4“  =  j  sin(M™S) 


th 


it 


cos(2  /ih  —  /imag)  sin(2/r/i  —  /jmag) 
2A  ^sin(2/i/i  —  nmag)  —  cos(2 fih  —  ^imag) 

sin(/imag  —  2/ih)  —  sin  fimag  cos (/imag  —  2/i/i)  —  cos  p.mag  \ 


77I0C  _ 

’  —  4A2  V  cos(/xmag  —  2 fih)  —  cos  /imag  —  sin(/imag  —  2 fih)  +  sin  /imag 


Local  error  indicators  for  the  modified  Airy  equation 


13.  Global  error  indicators  for  Magnus-4 


^global  =  5(0,  T)  -  5mag( 0,  T) 


n= 0 


n= 0 


=  n  (s™s(tn,tn+1) + Eir)  -  n  s™*(tn,tn+1) 


N-l 

N—l 


N-l 


=  ]T  Slmas(t„+1,  T)  £1“  Smag(0,  tn)  +  O  i  ||£locH2/ft 

n= 0 

Notice  that 

o 

Smag(0,T)  =  Xmag(t„)  exp(Dmag(t„))  (Xmag(t„)) 


-1 


n=Ar— 1 

=  X  exp(JDmag(0,  T))  X-1  +  h.o.t. 


therefore 


xglobal «  X 


■iV-l 


X  exP (Dmag( t„+i, T))  X-'Ep X exp (Dmag(0,  tn)) 


Ln=0 


X-1 


14.  Global  error  indicators  II 


For  example,  the  local  error  term  El2oc  contributes  globally 


h 2 

^global  _ 


rN-l 


sin  (/imag(4)) 


Ln=0 


0 

exp(-i<7n) 


exp(icrn)' 

0 


X-1  +  h.o.t. 


where 

(T  n  =  2t,jA  —  (T  +  /i)A  +  f  {tn)  + 

To  the  same  order  of  approximation,  this  corresponds  to  an  oscillatory 
integral  of  the  type 

T 

h  bdd(/i,  A)  f  g(t)  e~2lXt  dt  =  (9 

Jo 


Global  error  indicators  for  the  modified  Airy  equation 


15.  Precomputation 


where 


Kn  o  1=  (K0  +  AKi  +  /jK2)n  o  / 


n  k 


=  EE  E  Kaio...oKano/| 

k= 0  j= 0  ya&V(j,k,n) 
n  k 

=  E  E  A3 

fc=0  j=0 


aGV(j:k,n) 

V(j,  k,  n)  =  PermutationsjO, . . . ,  0, 1, . . . ,  1,  2, . . . ,  2} 

n—k  j  k—j 


Precomputation  for  the  Neumann  series 


N  n  k 

S%u(t;  A,ri  =  EEE  r^.»w  A5' 

n=0  fc=0  j=0 
TV  k 

=EEAiwA;"H 

&;=0  j=0 

where 

TV  AT 

AjSM  =  E  ri,M«  =  E  E  K»1  o  •  •  ■  o  K„„  o 

n=k  n=k  a£V(j,k,n) 

•  Similar  expansion  for  the  exponent  in  the  Magnus  series 

•  But:  Cannot  precompute  exponentiation 

•  Other  possible  splittings? 


16.  Application:  Reaction  diffusion  systems 


Ut  =  BUtt  +  cl\  +  F(U) 


Example:  Autocatalytic  two-component  system 


ut  =  8u££  +  cu^  —  uvm 

Vt  =  V££  +  CV(:  +  UVm 


Front-type  boundary  conditions 

(u,  v)  — >  (1,  0)  as  x  — »  — oo 
(w,  i>)  — >  (0, 1)  as  x  — >  +oo 

Traveling  wave  in  moving  frame 

u(t,  t)  =  uc(0 


17.  Stability  of  traveling  waves 


Perturbation  ansatz: 


u(i,t)  =  uc(i)  +  u(t)e>d 


Plugging  this  into  the  reaction-diffusion  system 

Ut  =  BU #  +  dJj:  +  F(U) 


yields 

XU  =[B%  +  cld(  +  DF{Uc{0 )]  U 
with  U(£)  — >  0  as  £  — >  ±oo. 


Reformulation 


where  Y  =  (U,  U^),  and 

A)  = 


Y'  =  A(.Z,\)Y, 

O  I 

B-^X-DFiU'd)))  -cB 


lm[A] 


Spectrum  of  the  linearized  traveling  wave  operator 


Complex  X  plane 


Re[XJ 


18.  The  Evans  function 


Limiting  systems 

A±{ A)  =  lim  A(£,  A) 

>±oo 

•  Assume  has  a  /.’dimensional  unstable  manifold 

•  A+  has  an  n  —  /.’-dimensional  stable  manifold 

•  Look  for  intersection  under  the  "evolution"  of  the  BVP 

Wronskian 

D(  A)  =  e“  ■  det  (Yf({;  A)  •  •  •  Yk~^  A)  Yt++1(£;  A)  •  •  •  Y+_kfo  A)) 

=  ■  (Yf  A  ■  ■  ■  A  Yk~)  A  (Y+ ,  A  ■  •  ■  A  Y+_,) 

=  e_  ('J  Tr-U;,:'A,d'  ■  (£/_(£;  A)  A  £/+(£;  A)) 

(Prefactor  ensures  ^-independence.) 


19.  Properties  of  the  Evans  function 

(Evans,  1975;  Alexander,  Gardner  &  Jones,  1990) 

•  Zeros  correspond  to  eigenvalues 

•  Order  of  the  zero  corresponds  to  algebraic  multiplicity 

•  Analytic  to  the  right  of  the  essential  spectrum 

•  Can  use  argument  principle  to  determine  number  of  zeros  in  the 
right  half  plane 

Numerical  issues 

•  May  want  to  rescale  solution  by  expected  exponential  growth 

•  Computing  second  most  unstable  eigenspace  is  numerically 
unstable 

— >  Projection  methods  or  Exterior  product  representation 


20.  Performance  comparison 

Near  the  ground  state 

•  Magnus  works  well 

•  Implicit/ Explicit  Runge-Kutta  perform  better 

In  the  stiff  regime  (A  =  i y  for  y  large) 

•  Magnus  works  well  in  the  interesting  regime 

•  Implicit  RK  performs  better 

•  Explicit  RK  cannot  be  used 

In  the  oscillatory  regime  (near  the  essential  spectrum) 

•  Magnus  is  the  best 

•  Explicit  RK  clearly  worse  than  implicit  RK 
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Error  behavior  in  the  oscillatory  regime 
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21.  Conclusions 


Magnus  integrators  offer 

•  Unconditional  stability 

•  Superior  performance  in  highly  oscillatory  regimes 

•  Possibility  of  a  priori  stepsize  control 

Would  I  recommend  Magnus? 

•  If  higher  bound  states  exist  or  are  suspected 

•  If  implementation  time  does  not  matter  (e.g.  for  library  codes) 

•  Generally,  if  good  splitting  methods  can  be  found  (— >  Jitse  Niesen) 

•  Discontinuous  potentials? 


